knitr::opts_knit$set(root.dir = 'C:/Users/ashua/OneDrive - VUMC/Work/05_OngoingProjects/Lifespan/KaidiReports/GAMLSS_HCP/')

library(gamlss)
library(splines)
library(table1)
library(magrittr)
library(plyr)
library(dplyr)
library(ggplot2)

1 Summary

In this project, we would like to implement the GAMLSS model to model the growth curves/charts of the brain using the HCP data.

For this specific report: The deviation score / centiles are calculated based the population-level curves. The variable Family is still included as random effects in the modeling to deal with the correlation between subjects, but when calculating the centiles for each subject, the population-level model (the one without random effects) is used.

2 Data

data = read.csv("HCP_PSYXSEC_DataHarmonized_220829.csv")
# some data cleaning
# data$Family[is.na(data$Family)] = paste0("temp", 1:sum(is.na(data$Family)))
# data <- data %>% group_by(Family) %>% slice_sample(n=1)
# nrow(temp) == temp$Family %>% unique %>% length

data$dx_group = data$Diagnosis
data$dx_group[data$Diagnosis != "Normal Control"] = "Schizophrenia Spectrum"
hcp = subset(data, Dataset == "HCP")

# Adding fake family ID to the subject from PSYXSEC
# ind = data$Dataset == "PSYXSEC"
# data$Family[ind] = paste0("f_", 1:sum(ind))
# # subset(data, Dataset == "PSYXSEC")

data$Family %<>% as.factor()

data$group = NA
data$group[data$Dataset == "HCP"] = "Normal - HCP"
data$group[data$dx_group == "Normal Control" & data$Dataset == "PSYXSEC"] = "Normal - PSYXSEC"
data$group[data$dx_group == "Schizophrenia Spectrum"] = "SZ Spectrum - PSYXSEC"

Note: The random intercepts/offsets for Family is not included. I noticed that after including the random offset Family, it becomes difficult for the models to converge. For most of the family units in the data, they only have one subjects participating in the studies. Moreover, there are 211 subjects have NA Family ID. This analysis ingored the possible correlation between individuals from the same family.

The random intercepts/offsets for Family are now included in the modeling of the population mean (\(\mu\)) and variance (\(\mu^2 \times \sigma^2\)). All healthy controls from the HCP studies have family ID, however, those from the PSYXSEC study (445 subjects) have missing family ID.

table(hcp$Family) %>% hist(main = "Num of subjects in each family")

# sum(is.na(data$Family[data$Dataset != "HCP"]))
table1(~ Age + Site + Race + eTIV + Sex + dx_group | Dataset, data = data)
HCP
(N=1952)
PSYXSEC
(N=445)
Overall
(N=2397)
Age
Mean (SD) 34.5 (19.9) 29.2 (10.5) 33.5 (18.7)
Median [Min, Max] 29.0 [5.60, 100] 25.1 [15.6, 60.5] 28.0 [5.60, 100]
Site
AMGH 135 (6.9%) 0 (0%) 135 (5.6%)
AUCLA 130 (6.7%) 0 (0%) 130 (5.4%)
AUMinn 163 (8.4%) 0 (0%) 163 (6.8%)
AWashU 174 (8.9%) 0 (0%) 174 (7.3%)
DHarvard 177 (9.1%) 0 (0%) 177 (7.4%)
DUCLA 102 (5.2%) 0 (0%) 102 (4.3%)
DUMinn 135 (6.9%) 0 (0%) 135 (5.6%)
DWashU 111 (5.7%) 0 (0%) 111 (4.6%)
YUWash 825 (42.3%) 0 (0%) 825 (34.4%)
Scan1 0 (0%) 280 (62.9%) 280 (11.7%)
Scan2 0 (0%) 2 (0.4%) 2 (0.1%)
Scan3 0 (0%) 163 (36.6%) 163 (6.8%)
Race
AA 272 (13.9%) 124 (27.9%) 396 (16.5%)
Other 316 (16.2%) 23 (5.2%) 339 (14.1%)
White 1364 (69.9%) 298 (67.0%) 1662 (69.3%)
eTIV
Mean (SD) 1510000 (208000) 1520000 (193000) 1510000 (206000)
Median [Min, Max] 1510000 [843000, 2190000] 1530000 [722000, 2080000] 1520000 [722000, 2190000]
Sex
F 1071 (54.9%) 148 (33.3%) 1219 (50.9%)
M 881 (45.1%) 297 (66.7%) 1178 (49.1%)
dx_group
Normal Control 1952 (100%) 211 (47.4%) 2163 (90.2%)
Schizophrenia Spectrum 0 (0%) 234 (52.6%) 234 (9.8%)

3 GAMLSS

Separate GAMLSS models are fitted to each outcome using the subjects from the HCP study, respectively. Generalized gamma (GG) distribution is specified for each of the outcome.

The final model for the trajectories of all measures are specified as below. We use R.VA as an example:

\[ \mathrm{R.VA} \sim GG(\mu, \sigma, \nu) \]

\[ log(\mu) = \beta_{\mu} + \beta_{\mu, age} \times ns(Age, df = 3) + \beta_{\mu, sex} \times Sex + \beta_{\mu, race} \times Race + \beta_{\mu, eTIV}\times I(neTIV/10000) + \alpha_{Family} \]

\[ log(\sigma) = \beta_{\sigma} + \beta_{\sigma, age} \times ns(Age, df = 3) + \beta_{\sigma, sex} \times Sex + \beta_{\sigma, race} \times Race + \beta_{\sigma, neTIV}\times I(neTIV/10000) + \alpha_{Family} \] \[ \nu = \alpha_{\nu} \]

# control function for gamlss
con_f <- gamlss.control(c.crit = 0.001, n.cyc = 300)

fit_gamlss <- function(y, full.data = data, con = con_f) {
  # select the variables needed. The function doesn't allow the existence of missing value in any variables even the variables are not included in the modeling.
  # Note: there is a bug in `predict.gamlss` -- when your data frame is called `data`, there will be a problem in predicting...
  
  full.data = full.data[, c(y, "Dataset", "dx_group", "Age", "Sex", "Race", "neTIV", "Family")] 
  
  # HCP subjects
  hcp = subset(full.data, Dataset == "HCP")

  num_obs = nrow(hcp)
  hcp = na.omit(hcp)
  num_obs_na.rm = nrow(hcp)
  cat("Notes:", num_obs - num_obs_na.rm, "observations were removed due to missingness. \n" )
  
  mu.form = as.formula(paste0(y, "~ ns(Age, df = 3) + Sex + Race + I(neTIV/10000) + re(random = ~1|Family)"))
  sigma.form = ~ ns(Age, df = 3) + Sex + Race + I(neTIV/10000) + re(random = ~1|Family)
  
  fit = gamlss(mu.form,
               sigma.formula = sigma.form,
               nu.formula = ~1,
               data = hcp,
               method = RS(), # or `CG()` or mixed()
               family = GG,
               control=con)
  return(fit)
}
# Test
#fit = fit_gamlss("nR.VA")
# @param centile_only: if only return centile estimates, no plots. Default = FALSE

plot_gamlss = function(obj, title, full.data = data, centile_only = FALSE){
  
  # 1. DATASET
  y = obj$mu.formula[[2]] %>% as.character()
  full.data = full.data[, c("Subject", y, "Dataset", "dx_group", "group", "Age", "Sex", "Race", "neTIV", "Family")] 
  
  # HCP subjects
  hcp = subset(full.data, Dataset == "HCP")
  # hcp = model.matrix(obj)
  
  # Schizophrenia Spectrum Subjects
  sz = subset(full.data, Dataset != "HCP")
  
  # 2. POPULATION CURVE (averaged over the sample)
  if (!centile_only){
    age_range = range(hcp$Age)
    age_grid = seq(from = age_range[1], to = age_range[2], by = 0.5)
    pred_data = hcp[, colnames(hcp) != "Age"]
    pred_data = cbind(Age = rep(age_grid, each = nrow(pred_data)), pred_data)
  
    # predict "mu" (without random effects)
    # Note: you have to add the original data here....
    mean_mu = predict(obj, data = hcp, newdata = pred_data, what = "mu", type = "response")
    # predict "sigma"
    mean_sigma = predict(obj, data = hcp, newdata = pred_data, what = "sigma", type = "response")
    # predict "nu"
    mean_nu = predict(obj, data = hcp, newdata = pred_data, what = "nu", type = "response")
    ## obtain the corresponding 2.5% and 97.5% centiles for mu for each subject
    mu_median =  gamlss.dist::qGG(0.5, mu = mean_mu, sigma = mean_sigma , nu = mean_nu )
    mu_LL = gamlss.dist::qGG(0.05, mu = mean_mu, sigma = mean_sigma , nu = mean_nu )
    mu_UL = gamlss.dist::qGG(0.95, mu = mean_mu, sigma = mean_sigma , nu = mean_nu )
    
    mean_var = mean_sigma^2 * mean_mu^2  
    
    pred_data = cbind(mean_mu, mu_median, mu_LL, mu_UL, mean_sigma, mean_var, mean_nu, pred_data)
  
    ## marginal means of median & variance at each age value
    marg_mean_mu_median = aggregate( mu_median ~ Age*Sex, data = pred_data, mean )
    marg_mean_mu_LL = aggregate( mu_LL ~ Age*Sex, data = pred_data, mean )
    marg_mean_mu_UL = aggregate( mu_UL ~ Age*Sex, data = pred_data, mean )
    marg_mean_sigma = aggregate(mean_sigma ~ Age*Sex, data = pred_data, mean )
    marg_mean_var = aggregate(mean_var ~ Age*Sex, data = pred_data, mean)
  
    # 3. PLOT
    xlim = c(0, 100)
    yrange = range(hcp[,y])
    ymin = round_any(yrange[1],10,f=floor)
    ymax = round_any(yrange[2],10,f=ceiling)
    ylim = c(ymin,ymax)
    ybreaks = floor((ylim[2]-ylim[1])/4)
    

    ## mu over age
    #main = paste0("Population trajectory of the median for ", y, "\n (with 5% and 95% centiles (dotted lines))"), 
    par(mar = c(4, 5.2, 2.5, 0.5)+0.03)
    #,mgp = c(2, 2, 0)) # Set the margin on all sides to 2
    
    colors <- c("pink", 
            "dodgerblue")

    plot(as.formula(paste0(y, "~ Age")), data = hcp, frame.plot = FALSE, 
         xlab = substitute(paste(bold("Age (years)"))),
         xlim = xlim, ylim = ylim,
         yaxt = "n",
         ylab = expression(bold(paste("Volumes (mm"^3*")"))), 
         main = title,
         col = colors[factor(Sex)],
         cex.lab=1.7,
         cex.main=2)
    
    #adjustcolor("dodgerblue", alpha.f = 0.35)) 
    
    axis(side = 2, at=seq(ylim[1],ylim[2],ybreaks))
  
    lines(mu_median ~ Age, data = marg_mean_mu_median[marg_mean_mu_median$Sex=="F",], type = "l", col = "deeppink", lwd = 4)
    lines(mu_LL ~ Age, data = marg_mean_mu_LL[marg_mean_mu_median$Sex=="F",], type = "l", lty = 3, col = "deeppink", lwd = 3)
    lines(mu_UL ~ Age, data = marg_mean_mu_UL[marg_mean_mu_median$Sex=="F",], type = "l", lty = 3, col = "deeppink", lwd = 3)

    lines(mu_median ~ Age, data = marg_mean_mu_median[marg_mean_mu_median$Sex=="M",], type = "l", col = "blue", lwd = 4)
    lines(mu_LL ~ Age, data = marg_mean_mu_LL[marg_mean_mu_median$Sex=="M",], type = "l", lty = 3, col = "blue", lwd = 3)
    lines(mu_UL ~ Age, data = marg_mean_mu_UL[marg_mean_mu_median$Sex=="M",], type = "l", lty = 3, col = "blue", lwd = 3)
 
    
    ## variance over age
    plot(mean_var ~ Age, data = marg_mean_var, type = "l", frame.plot = FALSE,
         lwd = 4, col = "dodgerblue",
         ylab = "Variance", xlab = "Age",
         main = paste0("Normative variance for ", y))
  }

    
  # 4. INDIVIDUAL DEVIANCE (i.e., deriving centile conditional on subject characteristics)
  # Note: function `centiles()` only works for model fits with one covariate.
  
  # 4.1 predict the centiles for HCP subjects 
  # ~~sing `predict` function to obtain subject-specific distribution parameter values~~
  # ~~Note: the predicted values are based on the values of the random effects.~~
  # NEW: the predicted values are NOT based on the random effects
  rdata = full.data
  ## design matrix for HCP subset
  ns_obj = ns(hcp$Age, df = 3)
  ns_age = predict(ns_obj, rdata$Age)
  sex = ifelse(rdata$Sex == "M", 1, 0)
  race_other = ifelse(rdata$Race == "Other", 1, 0)
  race_white = ifelse(rdata$Race == "White", 1, 0)
  neTIV_10000 = rdata$neTIV/10000
  mod_mat = cbind(1, ns_age, sex, race_other, race_white, neTIV_10000)
  
  ## coef estimates
  coef_all = coefAll(obj)
  mu_beta_hat = coef_all$mu[1:ncol(mod_mat)]
  sigma_beta_hat = coef_all$sigma[1:ncol(mod_mat)]
  nu_beta_hat = coef_all$nu
  ## predicted values
  pred_mu = exp(mod_mat %*% mu_beta_hat)
  pred_sigma = exp(mod_mat %*% sigma_beta_hat)
  pred_nu = rep(nu_beta_hat, nrow(rdata))
  rdata = cbind(rdata, pred_mu, pred_sigma, pred_nu)
  
  # pred_mu = predict(obj, data = hcp, what = "mu", type = "response")
  # pred_sigma = predict(obj, data = hcp, what = "sigma", type = "response")
  # pred_nu = predict(obj, data = hcp, what = "nu", type = "response")
  # hcp = cbind(hcp, pred_mu, pred_sigma, pred_nu)
  # obtain individual percentiles from subject-specific GG distributions
  rdata$subj_centile = gamlss.dist::pGG(rdata[, y], mu = rdata$pred_mu, sigma = rdata$pred_sigma, nu = rdata$pred_nu)
  
  ## DENSITY PLOT OF DEVIANCE
  
  if (!centile_only){
      plot_centile = ggplot(rdata, aes(x = subj_centile, fill = group, col = group)) + geom_histogram(aes(y=..density..), alpha = 0.2, bins = 30, position="identity") + geom_density(alpha=.2) + labs(title = paste("Centile density plot for", y), x = "Centiles", y = "Density") + facet_wrap(~group, ncol = 1)
    centile_med = aggregate( subj_centile ~ group, data = rdata, median)
    plot_centile = plot_centile + geom_point(data = centile_med, aes(x = subj_centile, y = 0, col = group), size = 3) + labs(title = y, x = "Centiles", y = "Density")
    print(plot_centile)

    # # individual deviance (centiles) vs rank
    temp_data = subset(rdata, Dataset != "HCP" & dx_group == "Normal Control")
    temp_data$rank_n = rank(temp_data$subj_centile)/nrow(temp_data)
    plot_centile_v_rank = ggplot(temp_data, aes(x = rank_n, y = subj_centile, col = group)) + geom_point() + geom_abline(slope = 1) + labs(title = paste("rank/sample size vs centiles for", y, "- NC from the PSYXSEC study"), x = "Rank scaled by sample size", y = "Centiles")
    print(plot_centile_v_rank)
  }
  # return centiles
  if (centile_only) return(rdata)
}
# TEST
# plot_gamlss(fit)

4 Model Results

4.1 AV

#fit_R.AV <- fit_gamlss(y = "nR.AV")
#saveRDS(fit_R.AV, file = "fit/fit_R.AV.rds")
fit_R.AV <- readRDS("gamlss_fit_obj/fit_R.AV.rds")
summary(fit_R.AV)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.336486   0.030901 140.334  < 2e-16 ***
## ns(Age, df = 3)1 -0.090390   0.017176  -5.262 1.61e-07 ***
## ns(Age, df = 3)2 -0.395426   0.030133 -13.123  < 2e-16 ***
## ns(Age, df = 3)3 -0.488810   0.023869 -20.479  < 2e-16 ***
## SexM             -0.005728   0.007590  -0.755   0.4505    
## RaceOther        -0.047234   0.011471  -4.118 4.02e-05 ***
## RaceWhite        -0.023870   0.009515  -2.509   0.0122 *  
## I(neTIV/10000)    0.004947   0.000198  24.989  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.704136   0.169256 -10.068  < 2e-16 ***
## ns(Age, df = 3)1  0.179870   0.084810   2.121  0.03408 *  
## ns(Age, df = 3)2  0.094036   0.170651   0.551  0.58168    
## ns(Age, df = 3)3  0.210751   0.118277   1.782  0.07496 .  
## SexM              0.116232   0.042607   2.728  0.00644 ** 
## RaceOther        -0.067651   0.059252  -1.142  0.25372    
## RaceWhite        -0.079765   0.048138  -1.657  0.09771 .  
## I(neTIV/10000)   -0.001974   0.001106  -1.785  0.07448 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.5945     0.3897   4.092 4.48e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  306.8979
##       Residual Deg. of Freedom:  1645.102 
##                       at cycle:  32 
##  
## Global Deviance:     16991.62 
##             AIC:     17605.41 
##             SBC:     19316.86 
## ******************************************************************
plot_gamlss(fit_R.AV,"Right Anteroventral")
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# con_temp <- gamlss.control(c.crit = 0.001, n.cyc = 300, mu.step = 0.1, sigma.step = 0.1, nu.step = 0.1)
# fit_L.AV <- fit_gamlss(y = "nL.AV", con = con_temp) 
# saveRDS(fit_L.AV, file = "fit/fit_L.AV.rds")
fit_L.AV <- readRDS("gamlss_fit_obj/fit_L.AV.rds")
summary(fit_L.AV)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.3178199  0.0321116 134.463  < 2e-16 ***
## ns(Age, df = 3)1 -0.0987395  0.0170980  -5.775 9.30e-09 ***
## ns(Age, df = 3)2 -0.3806025  0.0308155 -12.351  < 2e-16 ***
## ns(Age, df = 3)3 -0.4929875  0.0231753 -21.272  < 2e-16 ***
## SexM              0.0034004  0.0075974   0.448    0.655    
## RaceOther        -0.0621468  0.0115845  -5.365 9.35e-08 ***
## RaceWhite        -0.0557656  0.0096870  -5.757 1.03e-08 ***
## I(neTIV/10000)    0.0045443  0.0002009  22.623  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.324459   0.162321  -8.159 6.92e-16 ***
## ns(Age, df = 3)1  0.204368   0.083439   2.449   0.0144 *  
## ns(Age, df = 3)2 -0.051668   0.151582  -0.341   0.7333    
## ns(Age, df = 3)3  0.149183   0.108957   1.369   0.1711    
## SexM              0.182792   0.040628   4.499 7.34e-06 ***
## RaceOther        -0.079597   0.056698  -1.404   0.1606    
## RaceWhite        -0.048108   0.046974  -1.024   0.3059    
## I(neTIV/10000)   -0.004358   0.001049  -4.156 3.41e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.3380     0.3818   8.744   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  412.4777
##       Residual Deg. of Freedom:  1539.522 
##                       at cycle:  58 
##  
## Global Deviance:     16709.26 
##             AIC:     17534.22 
##             SBC:     19834.44 
## ******************************************************************
plot_gamlss(fit_L.AV,"Left Anteroventral")
## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded

## Warning in data.frame(..., check.names = FALSE): There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

4.2 VA

# fit_R.VA <- fit_gamlss(y = "nR.VA")
# saveRDS(fit_R.VA, file = "fit/fit_R.VA.rds")
fit_R.VA <- readRDS("gamlss_fit_obj/fit_R.VA.rds")
summary(fit_R.VA)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.1006038  0.0241937 210.824  < 2e-16 ***
## ns(Age, df = 3)1 -0.1236249  0.0134055  -9.222  < 2e-16 ***
## ns(Age, df = 3)2 -0.3393843  0.0225541 -15.048  < 2e-16 ***
## ns(Age, df = 3)3 -0.4436330  0.0187185 -23.700  < 2e-16 ***
## SexM             -0.0104138  0.0060917  -1.710 0.087558 .  
## RaceOther        -0.0283830  0.0084564  -3.356 0.000809 ***
## RaceWhite         0.0116368  0.0073222   1.589 0.112209    
## I(neTIV/10000)    0.0046895  0.0001607  29.187  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.1746573  0.1438867 -15.114  < 2e-16 ***
## ns(Age, df = 3)1  0.2481766  0.0854961   2.903  0.00375 ** 
## ns(Age, df = 3)2  0.3331504  0.1653229   2.015  0.04406 *  
## ns(Age, df = 3)3  0.1901254  0.1165958   1.631  0.10317    
## SexM             -0.0030207  0.0099475  -0.304  0.76142    
## RaceOther        -0.0971061  0.0493676  -1.967  0.04936 *  
## RaceWhite         0.0109048  0.0169633   0.643  0.52042    
## I(neTIV/10000)   -0.0013173  0.0007948  -1.657  0.09765 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -0.9922     0.4673  -2.123   0.0339 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  421.4835
##       Residual Deg. of Freedom:  1530.517 
##                       at cycle:  33 
##  
## Global Deviance:     19127.36 
##             AIC:     19970.32 
##             SBC:     22320.77 
## ******************************************************************
fit_R.VA %>% plot_gamlss("Right Ventral Anterior")
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

# fit_L.VA <- fit_gamlss(y = "nL.VA")
# saveRDS(fit_L.VA, file = "fit/fit_L.VA.rds")
fit_L.VA <- readRDS("gamlss_fit_obj/fit_L.VA.rds")
summary(fit_L.VA)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.0828254  0.0219352 231.721   <2e-16 ***
## ns(Age, df = 3)1 -0.1404496  0.0130618 -10.753   <2e-16 ***
## ns(Age, df = 3)2 -0.4085284  0.0215384 -18.967   <2e-16 ***
## ns(Age, df = 3)3 -0.5232676  0.0209615 -24.963   <2e-16 ***
## SexM              0.0007075  0.0064123   0.110   0.9122    
## RaceOther        -0.0197369  0.0086920  -2.271   0.0233 *  
## RaceWhite         0.0034756  0.0078843   0.441   0.6594    
## I(neTIV/10000)    0.0043196  0.0001602  26.964   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.3774124  0.1349913 -17.612  < 2e-16 ***
## ns(Age, df = 3)1  0.3204877  0.0857250   3.739 0.000192 ***
## ns(Age, df = 3)2  0.5223207  0.1587604   3.290 0.001026 ** 
## ns(Age, df = 3)3  0.5748047  0.1119727   5.133 3.23e-07 ***
## SexM              0.0170516  0.0308118   0.553 0.580067    
## RaceOther        -0.1488782  0.0594528  -2.504 0.012384 *  
## RaceWhite        -0.1164401  0.0478658  -2.433 0.015109 *  
## I(neTIV/10000)   -0.0002262  0.0007329  -0.309 0.757665    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.3881     0.5283  -4.521 6.67e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  494.9939
##       Residual Deg. of Freedom:  1457.006 
##                       at cycle:  38 
##  
## Global Deviance:     18419.8 
##             AIC:     19409.79 
##             SBC:     22170.18 
## ******************************************************************
fit_L.VA %>% plot_gamlss("Left Ventral Anterior")
## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded

## Warning in data.frame(..., check.names = FALSE): There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

4.3 VL

# fit_R.VL <- fit_gamlss(y = "nR.VL")
# saveRDS(fit_R.VL, file = "fit/fit_R.VL.rds")
fit_R.VL <- readRDS("gamlss_fit_obj/fit_R.VL.rds")
summary(fit_R.VL)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.431e+00  1.631e-02 394.423  < 2e-16 ***
## ns(Age, df = 3)1 -1.062e-01  1.130e-02  -9.397  < 2e-16 ***
## ns(Age, df = 3)2 -3.341e-01  1.695e-02 -19.709  < 2e-16 ***
## ns(Age, df = 3)3 -4.270e-01  1.521e-02 -28.063  < 2e-16 ***
## SexM             -3.651e-03  1.171e-03  -3.119  0.00185 ** 
## RaceOther        -1.442e-02  6.148e-03  -2.345  0.01919 *  
## RaceWhite         8.041e-03  5.643e-03   1.425  0.15439    
## I(neTIV/10000)    3.936e-03  9.319e-05  42.235  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.1289540  0.1243124 -17.126  < 2e-16 ***
## ns(Age, df = 3)1  0.1787333  0.1855895   0.963 0.335693    
## ns(Age, df = 3)2  0.5052426  0.2433606   2.076 0.038074 *  
## ns(Age, df = 3)3  0.4293012  0.1866583   2.300 0.021605 *  
## SexM              0.0109421  0.0117163   0.934 0.350513    
## RaceOther        -0.1572408  0.0682413  -2.304 0.021364 *  
## RaceWhite        -0.0956725  0.0350703  -2.728 0.006455 ** 
## I(neTIV/10000)   -0.0030741  0.0008991  -3.419 0.000648 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.29756    0.06088  -4.888 1.14e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  611.083
##       Residual Deg. of Freedom:  1340.917 
##                       at cycle:  31 
##  
## Global Deviance:     22859.79 
##             AIC:     24081.95 
##             SBC:     27489.73 
## ******************************************************************
fit_R.VL %>% plot_gamlss("Right Ventrolateral")
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

# fit_L.VL <- fit_gamlss(y = "nL.VL")
# saveRDS(fit_L.VL, file = "fit/fit_L.VL.rds")
fit_L.VL <- readRDS("gamlss_fit_obj/fit_L.VL.rds")
summary(fit_L.VL)
## Warning in summary.gamlss(fit_L.VL): summary: vcov has failed, option qr is used instead
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,  
##     family = GG, data = hcp, method = RS(), control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.4164555  0.0180096 356.280   <2e-16 ***
## ns(Age, df = 3)1 -0.1255590  0.0101464 -12.375   <2e-16 ***
## ns(Age, df = 3)2 -0.4018762  0.0166799 -24.093   <2e-16 ***
## ns(Age, df = 3)3 -0.4768878  0.0145035 -32.881   <2e-16 ***
## SexM             -0.0070763  0.0045983  -1.539   0.1240    
## RaceOther        -0.0167174  0.0066761  -2.504   0.0124 *  
## RaceWhite        -0.0038386  0.0057564  -0.667   0.5050    
## I(neTIV/10000)    0.0039653  0.0001202  32.993   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.168233   0.163836 -13.234  < 2e-16 ***
## ns(Age, df = 3)1  0.341323   0.085429   3.995 6.79e-05 ***
## ns(Age, df = 3)2  0.545283   0.161556   3.375 0.000758 ***
## ns(Age, df = 3)3  0.386741   0.115105   3.360 0.000800 ***
## SexM              0.039165   0.041883   0.935 0.349893    
## RaceOther        -0.132506   0.059284  -2.235 0.025566 *  
## RaceWhite        -0.087382   0.048322  -1.808 0.070767 .  
## I(neTIV/10000)   -0.003180   0.001075  -2.957 0.003153 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.5194     0.3969  -1.308    0.191
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms may not be reliable. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  540.4578
##       Residual Deg. of Freedom:  1411.542 
##                       at cycle:  30 
##  
## Global Deviance:     22688.86 
##             AIC:     23769.77 
##             SBC:     26783.7 
## ******************************************************************
fit_L.VL %>% plot_gamlss("Left Ventrolateral")
## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded

## Warning in data.frame(..., check.names = FALSE): There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

4.4 VPL

# con_temp <- gamlss.control(c.crit = 0.001, n.cyc = 300, mu.step = 0.01, sigma.step = 0.01, nu.step = 0.01)
# fit_R.VPL <- fit_gamlss(y = "nR.VPL", con = con_temp) 
# saveRDS(fit_R.VPL, file = "fit/fit_R.VPL.rds")
fit_R.VPL <- readRDS("gamlss_fit_obj/fit_R.VPL.rds")
summary(fit_R.VPL)
## Warning in summary.gamlss(fit_R.VPL): summary: vcov has failed, option qr is used instead
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,  
##     family = GG, data = hcp, method = RS(), control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.1381379  0.0215518 238.409  < 2e-16 ***
## ns(Age, df = 3)1  0.0024010  0.0116605   0.206   0.8369    
## ns(Age, df = 3)2 -0.0807364  0.0198881  -4.060 5.11e-05 ***
## ns(Age, df = 3)3 -0.1880015  0.0163009 -11.533  < 2e-16 ***
## SexM              0.0340769  0.0054827   6.215 6.25e-10 ***
## RaceOther         0.0477738  0.0080571   5.929 3.59e-09 ***
## RaceWhite         0.0150638  0.0068230   2.208   0.0274 *  
## I(neTIV/10000)    0.0032241  0.0001426  22.616  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.234192   0.163893 -13.632   <2e-16 ***
## ns(Age, df = 3)1  0.188653   0.085462   2.207   0.0274 *  
## ns(Age, df = 3)2  0.555743   0.161604   3.439   0.0006 ***
## ns(Age, df = 3)3  0.269718   0.115153   2.342   0.0193 *  
## SexM              0.040075   0.041898   0.956   0.3390    
## RaceOther        -0.094750   0.059306  -1.598   0.1103    
## RaceWhite        -0.071702   0.048342  -1.483   0.1382    
## I(neTIV/10000)   -0.001886   0.001076  -1.754   0.0797 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.01304    0.34109  -0.038    0.969
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms may not be reliable. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  500.9787
##       Residual Deg. of Freedom:  1451.021 
##                       at cycle:  36 
##  
## Global Deviance:     18503.05 
##             AIC:     19505.01 
##             SBC:     22298.77 
## ******************************************************************
fit_R.VPL %>% plot_gamlss("Right Ventral Posterolateral")
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

# 
# fit_L.VPL <- fit_gamlss(y = "nL.VPL")
# saveRDS(fit_L.VPL, file = "fit/fit_L.VPL.rds")
fit_L.VPL <- readRDS("gamlss_fit_obj/fit_L.VPL.rds")
summary(fit_L.VPL)
## Warning in summary.gamlss(fit_L.VPL): summary: vcov has failed, option qr is used instead
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,  
##     family = GG, data = hcp, method = RS(), control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.175952   0.022239 232.740  < 2e-16 ***
## ns(Age, df = 3)1  0.002724   0.012273   0.222   0.8244    
## ns(Age, df = 3)2 -0.112207   0.020710  -5.418 6.78e-08 ***
## ns(Age, df = 3)3 -0.220563   0.017225 -12.805  < 2e-16 ***
## SexM              0.026222   0.005690   4.608 4.33e-06 ***
## RaceOther         0.032999   0.008381   3.937 8.53e-05 ***
## RaceWhite         0.016827   0.007156   2.351   0.0188 *  
## I(neTIV/10000)    0.003234   0.000148  21.855  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.222516   0.163864 -13.563  < 2e-16 ***
## ns(Age, df = 3)1  0.199709   0.085445   2.337 0.019552 *  
## ns(Age, df = 3)2  0.540559   0.161578   3.346 0.000841 ***
## ns(Age, df = 3)3  0.298066   0.115129   2.589 0.009716 ** 
## SexM              0.051929   0.041890   1.240 0.215294    
## RaceOther        -0.117022   0.059294  -1.974 0.048608 *  
## RaceWhite        -0.081936   0.048332  -1.695 0.090221 .  
## I(neTIV/10000)   -0.001655   0.001076  -1.539 0.124115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.3206     0.3367  -0.952    0.341
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms may not be reliable. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  414.7198
##       Residual Deg. of Freedom:  1537.28 
##                       at cycle:  36 
##  
## Global Deviance:     18722.16 
##             AIC:     19551.6 
##             SBC:     21864.33 
## ******************************************************************
fit_L.VPL %>% plot_gamlss("Left Ventral Posterolateral")
## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded

## Warning in data.frame(..., check.names = FALSE): There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

4.5 Pul

# fit_R.Pul <- fit_gamlss(y = "nR.Pul")
# saveRDS(fit_R.Pul, file = "fit/fit_R.Pul.rds")
fit_R.Pul <- readRDS("gamlss_fit_obj/fit_R.Pul.rds")
summary(fit_R.Pul)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.9999426  0.0198245 353.095  < 2e-16 ***
## ns(Age, df = 3)1 -0.1540033  0.0113824 -13.530  < 2e-16 ***
## ns(Age, df = 3)2 -0.5157236  0.0184744 -27.916  < 2e-16 ***
## ns(Age, df = 3)3 -0.5236371  0.0165736 -31.595  < 2e-16 ***
## SexM              0.0357509  0.0051067   7.001  3.9e-12 ***
## RaceOther         0.0695703  0.0076612   9.081  < 2e-16 ***
## RaceWhite         0.0634289  0.0064723   9.800  < 2e-16 ***
## I(neTIV/10000)    0.0018047  0.0001318  13.689  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.7560376  0.1433865 -19.221  < 2e-16 ***
## ns(Age, df = 3)1  0.2710613  0.0835099   3.246 0.001198 ** 
## ns(Age, df = 3)2  0.6899750  0.1600834   4.310 1.74e-05 ***
## ns(Age, df = 3)3  0.4070112  0.1094739   3.718 0.000209 ***
## SexM              0.0073490  0.0250248   0.294 0.769055    
## RaceOther        -0.0989740  0.0592871  -1.669 0.095256 .  
## RaceWhite        -0.1359840  0.0475934  -2.857 0.004336 ** 
## I(neTIV/10000)    0.0012873  0.0008579   1.500 0.133715    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.3747     0.5348    6.31 3.72e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  519.4501
##       Residual Deg. of Freedom:  1432.55 
##                       at cycle:  31 
##  
## Global Deviance:     24211.35 
##             AIC:     25250.25 
##             SBC:     28147.02 
## ******************************************************************
fit_R.Pul %>% plot_gamlss("Right Pulvinar")
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

# fit_L.Pul <- fit_gamlss(y = "nL.Pul")
# saveRDS(fit_L.Pul, file = "fit/fit_L.Pul.rds")
fit_L.Pul <- readRDS("gamlss_fit_obj/fit_L.Pul.rds")
summary(fit_L.Pul)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.9976200  0.0203387 344.055  < 2e-16 ***
## ns(Age, df = 3)1 -0.1386545  0.0119464 -11.606  < 2e-16 ***
## ns(Age, df = 3)2 -0.4880280  0.0197342 -24.730  < 2e-16 ***
## ns(Age, df = 3)3 -0.4868830  0.0191067 -25.482  < 2e-16 ***
## SexM              0.0330419  0.0052252   6.324 3.34e-10 ***
## RaceOther         0.0749004  0.0078438   9.549  < 2e-16 ***
## RaceWhite         0.0773936  0.0065637  11.791  < 2e-16 ***
## I(neTIV/10000)    0.0020230  0.0001347  15.018  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.403251   0.162671 -14.774  < 2e-16 ***
## ns(Age, df = 3)1  0.238463   0.084560   2.820  0.00486 ** 
## ns(Age, df = 3)2  0.809308   0.159989   5.059 4.73e-07 ***
## ns(Age, df = 3)3  0.684229   0.113063   6.052 1.80e-09 ***
## SexM              0.041448   0.042030   0.986  0.32422    
## RaceOther        -0.066980   0.060658  -1.104  0.26967    
## RaceWhite        -0.110740   0.048557  -2.281  0.02271 *  
## I(neTIV/10000)   -0.001283   0.001082  -1.186  0.23571    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.7255     0.5318   5.126 3.34e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  412.4122
##       Residual Deg. of Freedom:  1539.588 
##                       at cycle:  17 
##  
## Global Deviance:     24444.56 
##             AIC:     25269.38 
##             SBC:     27569.24 
## ******************************************************************
fit_L.Pul %>% plot_gamlss("Left Pulvinar")
## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded

## Warning in data.frame(..., check.names = FALSE): There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

4.6 LGN

# fit_R.LGN <- fit_gamlss(y = "nR.LGN")
# saveRDS(fit_R.LGN, file = "fit/fit_R.LGN.rds")
fit_R.LGN <- readRDS("gamlss_fit_obj/fit_R.LGN.rds")
summary(fit_R.LGN)
## Warning in summary.gamlss(fit_R.LGN): summary: vcov has failed, option qr is used instead
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,  
##     family = GG, data = hcp, method = RS(), control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.2011170  0.0257284 163.287  < 2e-16 ***
## ns(Age, df = 3)1 -0.0555857  0.0146651  -3.790 0.000155 ***
## ns(Age, df = 3)2 -0.3102835  0.0249621 -12.430  < 2e-16 ***
## ns(Age, df = 3)3 -0.5386605  0.0232848 -23.134  < 2e-16 ***
## SexM              0.0178171  0.0066408   2.683 0.007359 ** 
## RaceOther         0.0015038  0.0094978   0.158 0.874211    
## RaceWhite        -0.0371695  0.0081972  -4.534 6.13e-06 ***
## I(neTIV/10000)    0.0025200  0.0001721  14.643  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.3151720  0.1638457 -14.130  < 2e-16 ***
## ns(Age, df = 3)1  0.1182908  0.0854323   1.385  0.16637    
## ns(Age, df = 3)2  0.6428284  0.1615588   3.979 7.25e-05 ***
## ns(Age, df = 3)3  0.5774035  0.1150985   5.017 5.87e-07 ***
## SexM             -0.0152082  0.0418851  -0.363  0.71659    
## RaceOther        -0.1596555  0.0592881  -2.693  0.00716 ** 
## RaceWhite        -0.0909612  0.0483254  -1.882  0.05999 .  
## I(neTIV/10000)    0.0003481  0.0010754   0.324  0.74619    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.3456     0.2830   1.221    0.222
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms may not be reliable. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  431.9813
##       Residual Deg. of Freedom:  1520.019 
##                       at cycle:  32 
##  
## Global Deviance:     14754.76 
##             AIC:     15618.72 
##             SBC:     18027.71 
## ******************************************************************
fit_R.LGN %>% plot_gamlss("Right Lateral Geniculate")
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

# fit_L.LGN <- fit_gamlss(y = "nL.LGN")
# saveRDS(fit_L.LGN, file = "fit/fit_L.LGN.rds")
fit_L.LGN <- readRDS("gamlss_fit_obj/fit_L.LGN.rds")
summary(fit_L.LGN)
## Warning in summary.gamlss(fit_L.LGN): summary: vcov has failed, option qr is used instead
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,  
##     family = GG, data = hcp, method = RS(), control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.1886188  0.0268914 155.761  < 2e-16 ***
## ns(Age, df = 3)1 -0.0706637  0.0154386  -4.577 5.01e-06 ***
## ns(Age, df = 3)2 -0.3468106  0.0254375 -13.634  < 2e-16 ***
## ns(Age, df = 3)3 -0.5853796  0.0222919 -26.260  < 2e-16 ***
## SexM              0.0369125  0.0069437   5.316 1.18e-07 ***
## RaceOther         0.0126113  0.0101088   1.248   0.2123    
## RaceWhite        -0.0138185  0.0083933  -1.646   0.0998 .  
## I(neTIV/10000)    0.0025617  0.0001792  14.292  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.238e+00  1.638e-01 -13.658  < 2e-16 ***
## ns(Age, df = 3)1  3.200e-01  8.542e-02   3.747 0.000186 ***
## ns(Age, df = 3)2  4.770e-01  1.615e-01   2.953 0.003196 ** 
## ns(Age, df = 3)3  3.497e-01  1.151e-01   3.038 0.002417 ** 
## SexM             -1.818e-02  4.188e-02  -0.434 0.664297    
## RaceOther        -3.015e-02  5.928e-02  -0.509 0.611107    
## RaceWhite        -5.679e-02  4.832e-02  -1.175 0.240041    
## I(neTIV/10000)   -2.055e-05  1.075e-03  -0.019 0.984758    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.4187     0.2778  -1.507    0.132
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms may not be reliable. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  363.5109
##       Residual Deg. of Freedom:  1588.489 
##                       at cycle:  20 
##  
## Global Deviance:     14953.17 
##             AIC:     15680.19 
##             SBC:     17707.35 
## ******************************************************************
fit_L.LGN %>% plot_gamlss("Left Lateral Geniculate")
## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded

## Warning in data.frame(..., check.names = FALSE): There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

4.7 MGN

# fit_R.MGN <- fit_gamlss(y = "nR.MGN")
# saveRDS(fit_R.MGN, file = "fit/fit_R.MGN.rds")
fit_R.MGN <- readRDS("gamlss_fit_obj/fit_R.MGN.rds")
summary(fit_R.MGN)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.8486749  0.0281790 136.580  < 2e-16 ***
## ns(Age, df = 3)1 -0.1012390  0.0157224  -6.439 1.56e-10 ***
## ns(Age, df = 3)2 -0.5185034  0.0271597 -19.091  < 2e-16 ***
## ns(Age, df = 3)3 -0.5793734  0.0232749 -24.893  < 2e-16 ***
## SexM              0.0300637  0.0072555   4.144 3.59e-05 ***
## RaceOther         0.0527740  0.0106932   4.935 8.79e-07 ***
## RaceWhite         0.0368115  0.0088847   4.143 3.59e-05 ***
## I(neTIV/10000)    0.0020486  0.0001869  10.961  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.2799009  0.1522172 -14.978  < 2e-16 ***
## ns(Age, df = 3)1  0.1152035  0.0816016   1.412 0.158199    
## ns(Age, df = 3)2  0.4628272  0.1615942   2.864 0.004233 ** 
## ns(Age, df = 3)3  0.3692912  0.1083389   3.409 0.000668 ***
## SexM             -0.0164405  0.0357480  -0.460 0.645648    
## RaceOther        -0.0524711  0.0581812  -0.902 0.367262    
## RaceWhite        -0.0822500  0.0474837  -1.732 0.083425 .  
## I(neTIV/10000)    0.0006995  0.0009383   0.745 0.456101    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9616     0.4276   4.587 4.82e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  262.0145
##       Residual Deg. of Freedom:  1689.986 
##                       at cycle:  26 
##  
## Global Deviance:     13187.95 
##             AIC:     13711.98 
##             SBC:     15173.13 
## ******************************************************************
fit_R.MGN  %>% plot_gamlss("Right Medial Geniculate")
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

# fit_L.MGN <- fit_gamlss(y = "nL.MGN")
# saveRDS(fit_L.MGN, file = "fit/fit_L.MGN.rds")
fit_L.MGN <- readRDS("gamlss_fit_obj/fit_L.MGN.rds")
summary(fit_L.MGN)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.8313748  0.0262356 146.037  < 2e-16 ***
## ns(Age, df = 3)1 -0.0993767  0.0146293  -6.793 1.55e-11 ***
## ns(Age, df = 3)2 -0.4568395  0.0248393 -18.392  < 2e-16 ***
## ns(Age, df = 3)3 -0.5206643  0.0202392 -25.726  < 2e-16 ***
## SexM              0.0251153  0.0067533   3.719 0.000207 ***
## RaceOther         0.0425273  0.0099682   4.266 2.10e-05 ***
## RaceWhite         0.0196793  0.0080555   2.443 0.014675 *  
## I(neTIV/10000)    0.0020520  0.0001746  11.756  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.4228414  0.1619408 -14.961  < 2e-16 ***
## ns(Age, df = 3)1  0.2084008  0.0837955   2.487  0.01298 *  
## ns(Age, df = 3)2  0.4420514  0.1620983   2.727  0.00646 ** 
## ns(Age, df = 3)3  0.2158756  0.1135599   1.901  0.05748 .  
## SexM             -0.0083584  0.0482355  -0.173  0.86245    
## RaceOther         0.0185452  0.0743993   0.249  0.80319    
## RaceWhite        -0.0457813  0.0558899  -0.819  0.41283    
## I(neTIV/10000)    0.0008301  0.0011769   0.705  0.48071    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.6791     0.4233   3.966 7.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  358.2205
##       Residual Deg. of Freedom:  1593.78 
##                       at cycle:  13 
##  
## Global Deviance:     12902.2 
##             AIC:     13618.65 
##             SBC:     15616.3 
## ******************************************************************
fit_L.MGN %>% plot_gamlss("Left Medial Geniculate")
## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded

## Warning in data.frame(..., check.names = FALSE): There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

4.8 CM

# fit_R.CM <- fit_gamlss(y = "nR.CM")
# saveRDS(fit_R.CM, file = "fit/fit_R.CM.rds")
fit_R.CM <- readRDS("gamlss_fit_obj/fit_R.CM.rds")
summary(fit_R.CM)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.4991495  0.0311324 144.517  < 2e-16 ***
## ns(Age, df = 3)1 -0.0859788  0.0179499  -4.790 1.82e-06 ***
## ns(Age, df = 3)2 -0.4876218  0.0299364 -16.289  < 2e-16 ***
## ns(Age, df = 3)3 -0.4529192  0.0255408 -17.733  < 2e-16 ***
## SexM              0.0391498  0.0081736   4.790 1.82e-06 ***
## RaceOther         0.0890777  0.0116438   7.650 3.37e-14 ***
## RaceWhite         0.0657958  0.0096864   6.793 1.53e-11 ***
## I(neTIV/10000)    0.0015034  0.0002095   7.177 1.07e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.531210   0.166052 -15.243  < 2e-16 ***
## ns(Age, df = 3)1  0.258450   0.082778   3.122  0.00183 ** 
## ns(Age, df = 3)2  0.492311   0.160946   3.059  0.00226 ** 
## ns(Age, df = 3)3  0.249255   0.108104   2.306  0.02125 *  
## SexM             -0.058628   0.042169  -1.390  0.16462    
## RaceOther        -0.040915   0.058948  -0.694  0.48772    
## RaceWhite        -0.046447   0.047547  -0.977  0.32878    
## I(neTIV/10000)    0.002861   0.001079   2.651  0.00811 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.2473     0.3835    5.86 5.57e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  288.2003
##       Residual Deg. of Freedom:  1663.8 
##                       at cycle:  24 
##  
## Global Deviance:     15969.16 
##             AIC:     16545.56 
##             SBC:     18152.74 
## ******************************************************************
fit_R.CM %>% plot_gamlss("Right Centromedian")
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

# fit_L.CM <- fit_gamlss(y = "nL.CM")
# saveRDS(fit_L.CM, file = "fit/fit_L.CM.rds")
fit_L.CM <- readRDS("gamlss_fit_obj/fit_L.CM.rds")
summary(fit_L.CM)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.5114604  0.0285677 157.922  < 2e-16 ***
## ns(Age, df = 3)1 -0.0729940  0.0165206  -4.418 1.06e-05 ***
## ns(Age, df = 3)2 -0.4496853  0.0268738 -16.733  < 2e-16 ***
## ns(Age, df = 3)3 -0.4193886  0.0228457 -18.357  < 2e-16 ***
## SexM              0.0362262  0.0074598   4.856 1.32e-06 ***
## RaceOther         0.0701914  0.0114136   6.150 9.84e-10 ***
## RaceWhite         0.0669973  0.0092582   7.237 7.20e-13 ***
## I(neTIV/10000)    0.0015892  0.0001913   8.308  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.555025   0.162582 -15.715  < 2e-16 ***
## ns(Age, df = 3)1  0.298060   0.082560   3.610 0.000316 ***
## ns(Age, df = 3)2  0.554606   0.162581   3.411 0.000663 ***
## ns(Age, df = 3)3  0.258942   0.111402   2.324 0.020233 *  
## SexM             -0.049854   0.040910  -1.219 0.223171    
## RaceOther         0.001783   0.011640   0.153 0.878302    
## RaceWhite        -0.118319   0.036501  -3.242 0.001214 ** 
## I(neTIV/10000)    0.002594   0.001074   2.415 0.015852 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9395     0.3757   5.163 2.74e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  394.9273
##       Residual Deg. of Freedom:  1557.073 
##                       at cycle:  26 
##  
## Global Deviance:     15821.92 
##             AIC:     16611.77 
##             SBC:     18814.13 
## ******************************************************************
fit_L.CM  %>% plot_gamlss("Left Centromedian")
## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded

## Warning in data.frame(..., check.names = FALSE): There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

4.9 MD

# fit_R.MD <- fit_gamlss(y = "nR.MD")
# saveRDS(fit_R.MD, file = "fit/fit_R.MD.rds")
fit_R.MD <- readRDS("gamlss_fit_obj/fit_R.MD.rds")
summary(fit_R.MD)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.1509257  0.0217291 283.074  < 2e-16 ***
## ns(Age, df = 3)1 -0.0534377  0.0130422  -4.097 4.40e-05 ***
## ns(Age, df = 3)2 -0.3788348  0.0218202 -17.362  < 2e-16 ***
## ns(Age, df = 3)3 -0.3209589  0.0195392 -16.426  < 2e-16 ***
## SexM             -0.0055777  0.0059875  -0.932    0.352    
## RaceOther         0.0625791  0.0083124   7.528 8.84e-14 ***
## RaceWhite         0.0500609  0.0068067   7.355 3.14e-13 ***
## I(neTIV/10000)    0.0015626  0.0001506  10.377  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -3.1913203  0.1482312 -21.529  < 2e-16 ***
## ns(Age, df = 3)1  0.1317723  0.0744351   1.770  0.07688 .  
## ns(Age, df = 3)2  0.8101773  0.1514237   5.350 1.01e-07 ***
## ns(Age, df = 3)3  0.4673210  0.0980506   4.766 2.06e-06 ***
## SexM             -0.0982921  0.0369703  -2.659  0.00793 ** 
## RaceOther         0.0054464  0.0291271   0.187  0.85169    
## RaceWhite        -0.0229580  0.0359510  -0.639  0.52319    
## I(neTIV/10000)    0.0042454  0.0009619   4.414 1.09e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.5742     0.8903   10.75   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  457.7676
##       Residual Deg. of Freedom:  1494.232 
##                       at cycle:  76 
##  
## Global Deviance:     21552.63 
##             AIC:     22468.16 
##             SBC:     25020.95 
## ******************************************************************
fit_R.MD  %>% plot_gamlss("Right Mediodorsal")
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

# fit_L.MD <- fit_gamlss(y = "nL.MD")
# saveRDS(fit_L.MD, file = "fit/fit_L.MD.rds")
fit_L.MD <- readRDS("gamlss_fit_obj/fit_L.MD.rds")
summary(fit_L.MD)
## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  gamlss(formula = mu.form, sigma.formula = sigma.form,  
##     nu.formula = ~1, family = GG, data = hcp, method = RS(),  
##     control = con) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.1434077  0.0277858 221.099  < 2e-16 ***
## ns(Age, df = 3)1 -0.0438666  0.0159427  -2.752    0.006 ** 
## ns(Age, df = 3)2 -0.4067019  0.0267378 -15.211  < 2e-16 ***
## ns(Age, df = 3)3 -0.3461348  0.0227581 -15.209  < 2e-16 ***
## SexM              0.0024296  0.0071370   0.340    0.734    
## RaceOther         0.0508062  0.0106681   4.762 2.09e-06 ***
## RaceWhite         0.0500736  0.0084953   5.894 4.60e-09 ***
## I(neTIV/10000)    0.0014357  0.0001861   7.714 2.17e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.558733   0.159370 -16.055  < 2e-16 ***
## ns(Age, df = 3)1  0.164944   0.083612   1.973  0.04870 *  
## ns(Age, df = 3)2  0.488547   0.162819   3.001  0.00274 ** 
## ns(Age, df = 3)3  0.308737   0.110398   2.797  0.00523 ** 
## SexM             -0.046958   0.040568  -1.158  0.24724    
## RaceOther         0.040273   0.058324   0.691  0.48997    
## RaceWhite        -0.056568   0.047855  -1.182  0.23736    
## I(neTIV/10000)    0.002360   0.001054   2.239  0.02532 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.942      0.432   6.811 1.38e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1952 
## Degrees of Freedom for the fit:  393.6556
##       Residual Deg. of Freedom:  1558.344 
##                       at cycle:  34 
##  
## Global Deviance:     21899.1 
##             AIC:     22686.41 
##             SBC:     24881.68 
## ******************************************************************
fit_L.MD  %>% plot_gamlss("Left Mediodorsal")
## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded

## Warning in data.frame(..., check.names = FALSE): There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction
## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 
## new prediction

5 Centile Mahalanobis Distance (CMD)

To create a composite score to indicate an individual’s overall deviation from the normative curve, CMD is calculated for each subject. CMD for the \(i\)-th subject is formalized as \[ D_i(x) = \sqrt{(x - \mu)^T S^{-1} (x-\mu) } \]

where \(x\) is the obtained centiles for each subject, \(\mu\) and \(S\) are the means and covariance of centiles among the normal controls (HCP subjects).

# adding predicted centiles with respect to each outcome to the data
outcome_var = paste0(rep(c("nR.", "nL."), times = 9), rep(c("AV", "VA", "VL", "VPL", "Pul", "LGN", "MGN", "CM", "MD"), each = 2))
output_data = data[, c("Subject", "Family", "Diagnosis", "dx_group", "Dataset", "group", "Age", "Sex", "Race", "Site", "neTIV", outcome_var)]

# obtaining the centile for each outcome
for (k in outcome_var){
  outcome = substring(k, 2)
  fit = get(paste0("fit_", outcome))
  centiles = plot_gamlss(fit, centile_only = TRUE)
  # subject centile for this outcome
  output_data = output_data[order(output_data$Subject), ]
  centiles = centiles[order(centiles$Subject), ]
  output_data = merge(output_data, centiles[, c("Subject", "subj_centile")], by = "Subject")
  colnames(output_data)[colnames(output_data) == "subj_centile"] = paste0("centile_", outcome)
}

# CMD
ind_row = output_data$group == "Normal - HCP"
ind_col = substring(colnames(output_data), 1, 8) == "centile_" 
S = cov(output_data[ind_row, ind_col])
mu = colMeans(output_data[ind_row, ind_col])

center_X = apply(output_data[, ind_col], 1, function(y) {y - mu}) %>% t()

CMD = apply(center_X, 1, function(y) sqrt(c(y) %*% solve(S) %*% c(y)))

output_data = cbind(output_data, CMD)
write.csv(output_data, "data_centiles_main.csv")
plot_CMD = ggplot(output_data, aes(x = CMD, fill = group, col = group)) + geom_histogram(aes(y=..density..), alpha = 0.2, bins = 50, position="identity") + geom_density(alpha=.2) + labs(x = "Centiles", y = "Density") 
# + facet_wrap(~group, ncol = 1)
CMD_med = aggregate( CMD ~ group, data = output_data, median)
plot_CMD = plot_CMD + geom_point(data = CMD_med, aes(x = CMD, y = 0, col = group), size = 3) + labs(title = "Density of CMD by groups", x = "Centiles", y = "Density")
print(plot_CMD)

plot_CMD = ggplot(output_data, aes(x = CMD, fill = dx_group, col = dx_group)) + geom_histogram(aes(y=..density..), alpha = 0.2, bins = 50, position="identity") + geom_density(alpha=.2) + labs(x = "Centiles", y = "Density") 
# + facet_wrap(~group, ncol = 1)
CMD_med = aggregate( CMD ~ dx_group, data = output_data, median)
plot_CMD = plot_CMD + geom_point(data = CMD_med, aes(x = CMD, y = 0, col = dx_group), size = 3) + labs(title = "Density of CMD by diagnosis", x = "Centiles", y = "Density")
print(plot_CMD)